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1. Introduction 

Relativistic shock fronts and currents sheets in relativistic flows play an important role 
in astrophysical models of gamma-ray bursts (for a review see Piran [1]), of jets and 
of pulsar winds (see for instance Michel [2] and Kirk [3]). The underlying plasma is 
probably composed of electrons, positrons and protons, whose temperature may be 
relativistic, i.e. comparable to their rest mass energy. 

We present an algorithm that computes the linear dispersion relation of waves in 
such plasmas within a Vlasov-Maxwell description. The method used is based on that 
presented by Daughton [I] for non relativistic Maxwellians, and involves explicit time 
integration of particle orbits along the unperturbed trajectories. We modify and extend 
this method by changing the manner in which the roots of the dispersion relation are 
located and adopt a fully relativistic approach, i.e. relativistic temperatures as well as 
relativistic drift speeds. 

Particular emphasis is given to the two-stream or Weibel instability [5]. This 
instability is very important in astrophysical processes because it is able to generate 
a magnetic field by pumping free energy from the anisotropic momentum distribution of 
an unmagnetised plasma or from the kinetic drift energy. There is an extensive literature 
about the Weibel instability. Although the dispersion curves has been found in some 
special cases such as, for example a water-bag distribution function, (Yoon [6]), or a fully 
relativistic bi-Maxwellian distribution function, (Yoon [7]), looking for an analytical 
expression for the dispersion relation for a given equilibrium distribution function is a 
complicated or even impossible task. The water-bag distribution is also the preferred 
profile to analyse magnetic field generation in fast ignitor scenarios, (Silva et al. [8]) or 
in relativistic shocks, (Wiersma and Achterberg [9j, Lyubarsky and Eichler [ID])- The 
Weibel instability in a magnetised electron-positron pair plasma has been investigated 
by Yang et al [11 J for the water-bag and for a smooth distribution function and a 
covariant description has been formulated by Melrose [12] and by Schlickeiser [13J. 
General conditions for the existence of the relativistic Weibel instability for arbitrary 
distribution functions are discussed in [TJ]. Wave propagation in counter-streaming 
magnetised nonrelativistic Maxwellian plasmas are studied in [151 IB]- The stability 
properties of a nonrelativistic Harris current sheet have been studied by Daughton p[] 
and also by Silin et al. [T7j. Streaming instabilities in relativistic magnetised plasmas 
and superluminous wave propagation are discussed by Buti [T8| [19] . 

In this paper, we investigate the relativistic Weibel instability for two counter- 
streaming Maxwellian distribution functions in an unmagnetised pair plasma. Our 
approach follows that of Zelenyi et al. [20] , who studied the tearing mode instabilities in a 
relativistic Harris current sheet, by integrating first order perturbations of the relativistic 
Maxwellian distribution function along unperturbed particle trajectories. Our algorithm 
is compared and checked against the dispersion curves for a single component magnetised 
plasma and for an unmagnetised plasma with counter-streaming components in the non- 
relativistic case. We present results for the growth rate of the Weibel instability in a 



3 



hot unmagnetised pair plasma consisting of counter-streaming relativistic Maxwellians. 

The paper is organised as follows. In Sec. [2j we present the full set of non- 
linear equations governing the motion in the pair plasma and describe the equilibrium 
of the two counter-streaming relativistic Maxwellians. Next, in Sec. [3j the full set 
of linearised Vlasov-Maxwell equations for the perturbed electromagnetic potentials 
around the equilibrium state are presented. The eigenvalue problem and the algorithm 
are discussed in Sec. HJ Results for magnetised plasma wave oscillations and for the 
non-relativistic Weibel instability as well as for the relativistic Weibel instability are 
shown in Sec. |5j Conclusions are drawn in Sec. [61 



2. Equilibrium 

Our purpose is to study the Weibel instability in a relativistic plasma. This plasma is 
made of counter-streaming electrons and positrons with relativistic temperatures and 
evolving in a static external magnetic field aligned with the z-axis such that 

B = B {x) e z (1) 

In equilibrium, there is no electric field, Eq = and the charges drift in the y direction 
at a relativistic velocity ±U S , (+ for positrons and — for electrons with U s > 0). We use 
Cartesian coordinates, denoted by (x, y, z), and the corresponding basis (e x , e y , e z ). The 
distribution function at equilibrium for each species " s" , denoted by f 0s (r , p) , is assumed 
to be a relativistic Maxwellian with drift speed ±U S . Adopting the usual notations, 
namely t, r, v, p, m s , q s for respectively the time, position, 3- velocity, 3-momentum, mass 
and charge of a particle of species s, the stationary distribution function reads : 

fos(r,p) = 3 3 W n ( % M /O S eXp [~ r - & - UsP y )/&smsC 2 } (2) 

no s {r) is the particle number density, E the total energy of a particle, p y the y-component 



of its momentum, c the speed of light, Y S = \/\J\ — U^/c 2 the Lorentz factor associated 
with the drift motion and K 2 the modified Bessel function of the second kind. The 
temperature of the gas is normalised to the rest mass energy of the leptons such that 

e s = (3) 

and kg is the Boltzmann constant. We introduce the standard electromagnetic scalar 
and vector potentials (<j>,A), related to the electromagnetic field {E,B) by : 

^-W-^ (4) 

B = VAl (5) 
We employ the Lorenz gauge condition by imposing 

dwA + eofio^ = (6) 
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with e p c = 1. The relation between potentials and sources then reads : 

The source terms represented by the charge p and current j densities, are expressed in 
terms of the distribution functions by : 



J (r ; t )=Y^q s jjj -^-J^ p, t) d V (8 6) 

where 7 = -^/l + p 2 /m 2 c 2 is the Lorentz factor of a particle. The time evolution of 
the distribution functions f s is governed by the well-known relativistic Vlasov-Maxwell 
equations written for each species : 

f + -§ + 9 .(^AB).§| = (9) 

The self-consistent non-linear evolution of the plasma is entirely determined by the set 
of equations 

3. Linearisation 

Our next step is to investigate the stability properties of the equilibrium configuration 
given in (jSJ) • The Vlasov equation ([9]) is linearised for each species about the equilibrium 
state fos, ©, to obtain the time evolution of the first order perturbation as 

t = -?^ + ^^'| ( 10 ) 
Perturbations are denoted by the subscript 1 whereas the equilibrium quantities are 
depicted by a subscript 0. The total time derivative d/dt is to be taken along the 
trajectories of the particles in the unperturbed electromagnetic field (E = 0, B ). More 
explicitly, performing the time integration in ( flOl) . the perturbed distribution function 
is given by 

f ls (r,p,t) = -q f [^(rV) +v' A B^r ',*')] ~if,f)^ (11) 

We assume that the perturbation of the distribution function at initial time vanishes, i.e., 
fis{r,p, t = —00) = 0. Position and momentum are determined by solving the equations 
of motion for a single particle in the unperturbed electromagnetic field (E = 0,B ). 
Thus, the unperturbed trajectories are solutions of the following system of ordinary 
differential equations 

% - v'(t') (12«) 
^- = q v'(t')AB (r>,t>) (12b) 
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with initial conditions r'(t' = t) = f and p'ii! = t) = p. Note also that for the relativistic 
Maxwellian, (jSj), we have 

dfo.s T s fos 



(v-U s e Y ) 



(13) 



Expanding all physical scalar quantities if) in terms of plane waves with complex 
frequency uj and real wavenumber vector k contained in the plane (yOz) 

4>(r,t) = i/j(x) exp(i (k y y + k z z - tut)) (14) 

we arrive by standard techniques at expressions for the charge density and current 
density 



p(?,t)= 



6 s c 2 



[U a A y (r, t) - <p(r, t)+i(u-ky U s ) x 



x 




r(t) 



fo.fr P) / [—■A(f',r') 

m, 



V' 
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r ,r 



dr d p 



J(r, t) 



e s c 2 



+ i(u — ky U s ) 



(U s A y (r,t)-(f>(r,t)) 




pc 



(15 a) 
f 0s (r,p)d 3 p+ 




pc 



E 



fos(r,p) x 



r(t) 



— • A(r r') 

m s 



7 



f t') ) dr d 3 p 



;i56) 



where the Lorentz factor of an unperturbed trajectory is 7' = \/l + p' 2 /m 2 s c 2 . 
Liouville's theorem 

fos(f'(T'),p'(T'),r') = f 0s (f,^ (16) 

stating that /o s is constant along the unperturbed trajectories defined below by fll7oj) and 
( 11761) . was used to extract fo s from the final integration over r' . To compute the integrals, 
it is more convenient to use the proper time defined by dr' = dt'/'y'. This relation 
between proper and observer time also defines the limit of time integration r(t). The 
(non- relativistic) plasma frequency corresponding to species "s" is u 2 = n 0s q 2 /m s e 
and the normalised distribution function is defined by fo s = n 0s f 0s . The trajectories 
are integrated over the unperturbed orbits, in the proper frame, following the equations 
of motion flT^j) . (fl26j) : 

dr" _ p'{r') 
dr' m s 

with initial conditions r'ir' = r) = r and p'ir' = r) = p. The non-relativistic cyclotron 
frequency is given by ub s = qsB Q /m s . According to equation fl!7aj) and (I176j) . both 



(17a) 
(176) 
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energy 7 ' and momentum component p' z are conserved so that the equations of motion 
in the z-direction can be integrated analytically: 

z'{t') = z +p z r'/m s (18a) 

p' z {r') = p z = constant (18 6) 

Consequently, the p z integration in (I15al) and (11561) can be factored out and performed 
separately After some algebraic manipulations, the charge density is expressed as 

n 0s q 2 s T s 



x 



k B T s 

F 0s {x, p x , P y ) S^(x, p x , p y ) dp x dp y 



R 2 



and the current density as 



7*0*0 = E 



X 



^0*0 = Yl 



k B T s 
n 0s q 2 s r 



i (a; — A; y f/ s ) x 



^os {x,Px, p y ) Sj (2, p x , p y ) dp x dp y 



[T S U 2 Ay(x)+i(U-kyU S ) X 



X 



— F 0s (a;, p x , p y ) SjOc, p x , p y ) dp x dp y 
n 0s q 2 s T s . 



i (a; — k y U s ) x 



x 



fl 2 



F 0s (x,p x ,Py) S jz (x,p x ,p y ) dpxdpy 



For brevity, we introduced the following functions : 

exp(r s [/ s p y /6 s m s c 2 ) 



F 0s (a;,p x ,p y ) 

5^(x,p x ,Py) 
Sj[X, p x , Py) 
S jx (x,Px,Py) 



4vrm3 c 3 9^(1/0 



(pUx+P y ^)^ L + ^ J ( 



<//, 



Pi 



(pU;+p;4)%i + ^/,. 49 



(p',K+PyA' y )^ + A' z I 3 



(19a) 



(196) 



(19c) 



m 6 



(19d) 

(20 a) 

e^dr' 

(206) 
e^dr' 
(20c) 
^ dr' 
(20 d) 



We employ a short-hand notation for (0', A4, Ay, A^.) which should be understood as the 
scalar and vector potential evaluated at (r ', r') on the unperturbed trajectories ( 117ol) 

and (11761) . We also use the following relation for the relativistic Maxwellian (jSJ) 

-,2 




pc 



f 0s d p = T s c(3 s 



(21) 
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It is easy to show that the proper time integration along the unperturbed orbits in the 
z direction, (\18a\) and (j 1 8 &[) . leads to expressions involving modified Bessel functions of 
the second kind Kq, K\ and K 2 such that 



m s cy± 



m s cy ± 




ft 



a 



•■J A, 



m s c 7 ± 



3 2 
m s c 7 ± 




K 2 {2 ^Jaft) -K (2 y/ap) 



(22 a) 

(22b) 

(22c) 
(22 d) 

(22 e) 
(22/) 



a 

1 

The argument of the modified Bessel functions for wave propagation along e z (k = k z e z ) 
are given by : 

(23) 



P± 

7.1 
a 



\Jvl+p- 



1 + p^/m 2 c 2 



XL 
2 

7-L 



A/r s /@ s + z (u; - fc z c) r' 



T s /Q s + i (uj + k z c) 



(24) 
(25) 
(26) 



The integrals (l22al) -( 22^ ) were first derived by Trubnikov [21]. Note also that 
the above integrals are computed assuming a relativistic Maxwellian distribution 
function. The analytical integration over p z leads to expressions involving modified 
Bessel functions K n . In the general case of non-Maxwellian distribution functions, an 
analytical integration over p z might not be possible, in which case it would have to be 
performed numerically as for the p Xjy components. This will of course decrease the speed 
of computation of the dispersion relation. 



4. The eigenvalue system 

4-1. Derivation 

The eigenvalue system is found by solving the equations for the electromagnetic potential 
determined according to the source distribution given by (119 a)) , (1196j) . (11 9 eft and (119d|) . 
Inserting the latter expressions into (!7aj) and (ITfej) . the eigenvalue system reads : 



J"(aO - k 



A"(x) 



UJ 



UJ' 



<j>(x) 
A{x) 







(27a) 
(276) 
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where prime ' denotes derivative with respect to x. The set of equations (I27al) and (12761) 
can be written in terms of the unknown 4-dimensional vector \I/ = (<p, A) 

M(uj, k) • # = (28) 

It is a non-linear eigenvalue problem for the matrix M with eigenvector \1/ and 
eigenvalue u. Daughton [4J solved this system by locating the zeroes of the determinant 
of M. We modify his method slightly by solving simultaneously for both the eigenvalues 
and the eigenvectors. A Newton- Raphson algorithm is used to find to and To check if 
the matrix M is singular, or in other words, if the iteration process has converged, instead 
of computing the determinant of M, which should vanish, we attempt to zero the ratio 
between the smallest and largest of the eigenvalues Aj of the matrix linear eigenvalue 
equation M ■ u — Xu. This latter eigenvalue equation should not be confused with 
our original initial physical eigenvalue problem (1281) and is absolutely not related to the 
computation of the dispersion relation. This procedure of computing the eigenvalues A 
of M is simply another mean to verify if M is singular. It helps to decide whether the 
algorithm has converged or not. More details are given in the next subsection. This 
procedure allows us to track the dispersion curves through crossing points. An example 
is given below. 

4-2. Algorithm 

Finding the roots of the determinant of the matrix M is the traditional way to search 
for the non-linear eigenvalues in (12 8 j) . However, it is not the method of choice in our 
problem because of the existence of different branches of the dispersion relation. For 
multiple branches in the dispersion relation, it is more efficient to look simultaneously 
for the eigenvalues and the eigenvectors. 

As an example, consider a homogeneous plasma (extension to an inhomogeneous 
configuration is straightforward). In this case, the matrix M is of dimension 4x4, and 
\l/ has 4 components, one for <ft and three for A, each of which is complex. 

At this stage there are 10 unknowns: 

• the complex eigenvalue u (real and imaginary parts) ; 

• the complex eigenfunction \l/ (real and imaginary parts each of the 4 components). 

We first reduce the system by normalising the eigenvector to unity, 1 1 \& 1 1 = 1 and by 
setting its phase such that the imaginary part of the last component of \I> vanishes, 
which we refer to as phase locking. The remaining 8 independent unknowns must then 
be found by solving the nonlinear system (1281) . that consist of 8 equations. To do this, 
we use the standard, globally convergent Newton- Raphson method, [22J. 

For fixed k, the iteration starts with an initial guess for the eigenvalue u = uj and 
the eigenvector \I> = \I>o- An improved guess is found by Newton- Raphson iteration in 
8 dimensions. Each step involves evaluating the Jacobian of the system (1281) and uses 
line searches and backtracking as described in [22]. The iteration is stopped when the 
matrix M is close to being singular, i.e. when detM ~ 0. This condition is verified 



9 



indirectly by computing the four complex eigenvalues of M, (Ai, A 2 , A 3 , A 4 ), that are 

solutions of the linear eigenvalue problem M -u— Xu. We emphasize that the Aj should 

not be confused with the eigenvalues u of our original nonlinear problem. When the 

ratio of the smallest to the largest Aj is sufficiently small, namely when 

mim e [i 41 A; . . 

el j ' < e 29 

maxi e [i..4] ||Aj|| 

with, typically, e = 10~ 6 ...10~ 8 , we assume that convergence has been achieved and 
terminate the iteration procedure. 



5. Results 

We apply our algorithm to the following cases : 

• the familiar electromagnetic oscillation modes of a nonrelativistic magnetised 
plasma ; 

• a counter-streaming pair-plasma with nonrelativistic temperature and nonrelativis- 
tic drift speed corresponding to two oppositely shifted Maxwellian distribution func- 
tions ; 

• a counter-streaming pair plasma with relativistic temperature and drift speeds close 
to the velocity of light, U s = 0.05 — 0.3 c. 



5.1. Plasma oscillations 

The dispersion relation for a single species plasma in a homogeneous magnetic field 
obtained using a nonrelativistic version of our code is shown in figure [TJ The plasma 
frequency is normalised to unity u v — 1 and the thermal speed of the non-relativistic 
Maxwellian is v t h = 0.1 c. In figured} we compare our numerical results with the analytic 
expressions found using a simple root finding algorithm from the exact dispersion 
relation for wave propagation parallel to the magnetic field Bo, (k y = 0) in order to 
check the correctness of our algorithm implementation. For electrostatic waves, we 
found the usual dispersion relation 

.2 r / \ i 

(30) 



1 + 2 1 



where the plasma dispersion function Z is defined for Im(£) > by 

1 /" + °° e~* 2 

V 71 J-oo t - C 



and analytically continued for Im(£) < 0, see for instance Delcroix and Bers [23]. For 
the left and right-handed circularly polarised electromagnetic wave we found 

^ l + -^-z(^) (32) 



u 2 ujk z v th V k * v th 

The numerical results are compared with the analytical exact expression for 
electrostatic waves (green curve, (}3T)j) ) and electromagnetic waves (red and blue curves, 
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Plasma oscillations 
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Figure 1. Dispersion relation for the transverse and longitudinal modes in a non- 
relativistic magnetised plasma. The exact analytical expressions for electrostatic waves 
are shown in green curve and those for electromagnetic waves are shown in red and 
blue curves. The plasma frequency is normalised to unity. 



( 132]) ). The advantages of implementing our modified root finding algorithm by looking 
for eigenvalues simultaneously with eigenfunctions is evident. The algorithm has no 
difficulty tracking only one branch of the dispersion relation even when crossing another 
branch. However it requires more computational effort since the root finding occurs in 
a space of higher dimension than is needed to find the eigenvalue alone. 

To test our algorithm in cases where the imaginary part lm(u) is small compared to 
the real part Re(u;), we consider electrostatic plasma oscillations, in the short wavelength 
limit, k <C 1. Numerical results for (Re/im)(u;), respectively red and blue dots, are 
compared with the analytical dispersion relation, solid black curves, (1301 . figure [21 As 
long as the ratio Im(a;)/Re(a;) is larger than the accuracy, which is about 5 digits, the 
results are reliable. 



5.2. "Classical" Weibel instability 



Next we check our algorithm for a situation with two counter-streaming species of non- 
relativistic temperature, taking v^/c = 10~ 3 and a small drift speed, U s = 0.1c. 

We recall that the exact analytical dispersion relation for the classical counter- 
streaming Weibel instability for equal and opposite drift speed for both species with the 
same plasma frequency u> p is given by 



kU 2 



2 r 



1 



1 



ut 



1 + 



to 



(33) 



The results are adapted from Delcroix and Bers [23] • In figure [3] we show the dispersion 
relation for the Weibel instability (red points) compared with the analytical expression 
(blue curve, (1331 ) and the asymptotic matching in case of a zero temperature plasma 
(green lines). The agreement is excellent. 
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Plasma oscillations 




Figure 2. Dispersion relation for electrostatic plasma waves. The real and imaginary 
parts of the cigenfrequency are shown in red and blue dots respectively. The plasma 
frequency is normalised to unity. 
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Figure 3. Growth rate for the non-relativistic Weibel instability (<d s -C 1). Results 
from our algorithm are depicted with red points, the exact analytical expression in blue 
curve and the asymptotic matching in case of a zero temperature plasma (0 S = 0) in 
green lines. 



Results for the dispersion relation in the case of small drift speeds are also shown 
in figure HI Here, the large k z asymptote is absent due to thermal effects. 

5.3. Relativistic Weibel instability 

Finally we computed the dispersion relation for a two-component fully relativistic 
counter-streaming plasma with high temperature, G s = 10 2 and different drift speeds 
as shown in figure El The dispersion curves are similar to those obtained in the classical 
Weibel instability. Due to the spread in momentum present at finite temperature, the 
instability is suppressed for large wavenumbers k z in both the classical and relativistic 
cases. 
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Weibel instability 
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Figure 4. Growth rate for the non-relativistic Weibel instability (9 S <C 1) for different 
drift speeds U s . 
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Figure 5. Growth rate for the relativistic Weibel instability (0 S 3> 1) for different 
relativistic drift speeds U s . 

6. Conclusion 

We have constructed an algorithm to solve the dispersion relation for non-relativistic 
and relativistic multi-component plasmas. The code has been validated by comparing 
the results with some typical configurations with homogeneous magnetic field, for which 
the analytical dispersion relation is known. New results have been obtained for the 
growth rate of the relativistic Weibel or two-stream instability for two counter-streaming 
relativistic Maxwellians which complement those found using a water-bag distribution. 
This code is easily extended to inhomogeneous plasmas, and is therefore a suitable tool 
for the study of stability properties of configurations of interest in gamma-ray burst and 
pulsar wind theories. 
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